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Abstract 

The composition of the lower mantle can be investigated by examining densities and seismic 
velocities of compositional models as functions of depth. In order to do this it is necessary to 
know the volumes and thermoelastic properties of the compositional constituents under lower 
mantle conditions. We determined the thermal equation of state (EOS) of MgSiOs perovskite 
using the nonempirical variational induced breathing (VIB) interatomic potential with molecular 
dynamics simulations at pressures and temperatures of the lower mantle. We fit our 
pressure-volume-temperature results to a thermal EOS of the form P(y,T) = ^0(^5^0) + 
APth{T), where Tq = 300 K and Pq is the isothermal Universal EOS. The thermal pressure APth 
can be represented by a linear relationship APth = a + bT. We find Vq = 165.40 A^, K^^ = 273 
GPa, Klp^ = 3.86, a = -1.99 GPa, and b = 0.00664 GPa K"^ for pressures of 0-140 GPa and 
temperatures of 300-3000 K. By fixing Vq to the experimentally determined value of 162.49 
and calculating density and bulk sound velocity profiles along a lower mantle geotherm we find 
that the lower mantle cannot consist solely of (Mg,Fe)Si03 perovskite with Xng ranging from 
0.9-1.0. Using pyrolitic compositions of 67 vol % perovskite (XMg = 0.93-0.96) and 33 vol % 
magnesiowiistite (XMg = 0.82-0.86), however, we obtained density and velocity profiles that are 
in excellent agreement with seismological models for a reasonable geotherm. 
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1. Introduction 

By comparing density and seismic velocity profiles 
of compositional models with seismological models 
it is possible to investigate the composition of the 
Earth's lower mantle. Therefore an understanding of 
the high pressure, high temperature properties, and 
equation of state (EOS) of (Mg,Fe)Si03 perovskite is 
vital to our understanding of the lower mantle, as this 
mineral accounts for perhaps two thirds of the miner- 



alogy of this region [Bina, 1998] and one third of the 
volume of the entire planet. Experiments have been 
performed up to pressures approaching those found 
at the base of the mantle, but direct coverage of the 
lower mantle geotherm has been limited to perhaps 
the uppermost one third (Figure |^) . Experiments in- 
clude resistively heated mult ianvil apparati |iyawff ei 
al, 1994|; ^ato et al l"995| ; \Utsumi et al, 1995| ; \Fu- 
namori et al., 1996|, diamond anvil cells (DACs) at 
ambient temperature \Kudoh et al., 1987; Ross ana 



Hazen, 1990 1, resistively heated DACs [Mao et al 



1991; Saxena et al., 1999], laser heated DACs [Knittle 



and Jeanloz 1987; Piquet et al., 1998, 2000 1, and shock 



and Ahrens, 1999| 



wave experiments [Duffy and Ahrens, 1993; Akin. 



Theoretical methods, such as those used here, per- 
mit the investigation of thermodynamic and ther- 
moelastic properties of minerals at high pressures 
and temperatures. Previous theoretical work on or- 
thorhombic MgSiOa perovskite (Mg-pv) has been 
done using molecular dynamics (MD) using the semi- 
empirical potential of Matsui [1988] at 0-1000 K and 
0-10 GPa (at 300 K) iMatsui, 198^ and 300-5500 K 
and 0-400 GPa [Belonoshko, 1994 ; Belonoshko ana 
Dubrovinsky, 199t[ and combined with lattice dynam- 
ics at 500-3000 K and 0-100 GPa p'atel et al, 1996 ]. 



Wolf and Bukowinski [1987] used a rigid ion modi- 
fied electron gas (MEG) model at 0-800 K and 0-150 
GPa, while Hemley et al. J1987] used a MEG model 
combined with shell stabilized ion charge densities at 
0-2500 K and 0-200 GPa. Cohen [1987] did quasi- 
harmonic lattice dynamics calculations from 0-3000 K 
and 0-150 GPa using the potential induced breathing 
(PIB) model. First principles static lattice calcula- 
tions (T = K) have also been done up to pressures of 
150 GPa, using plane wave pseudopotential (PWPP) 
MD Wentzcovitch et al., 1995| ] and the linearized 
augmented plane wave (LAPW) method [Stixrude 



and Cohen, 1993 ]. Karki et al. [1997] also used the 
PWPP method to examine the athermal elastic mod- 
uli of Mg-pv. 



In addition to the equation of state the thermody- 
namic stability of perovskite in the lower mantle is an 
open question and has been studied with experimen- 
tal and theoretical methods. Work by Meade et al. 
[1995], Saxena et al. [1996, 1998], and Dubrovinsky et 
al. [1998] indicate that perovskite will break down to 
oxides under lower mantle conditions. The thermo- 
dynamic analysis of Stixrude and Bukowinski [1992] 
and the experimental work of Serghiou et al. [1998], 
however, suggest that it will not. 

We used MD simulations using the nonempirical 
variational induced breathing (VIB) model, similar to 
the model of Wolf and Bukowinski [1988], to investi- 
gate the properties and EOS of Mg-pv at pressures (0- 
140 GPa) and temperatures (300-3000 K) that cover 
the bulk of the conditions found in the lower mantle. 
Newton's equations of motion are solved as functions 
of time, and equilibrium properties are obtained as 
time averages over a sufficiently long interval. This 
accounts for all orders of anharmonicity but not for 
quantum lattice dynamics effects. Thus our results 
are more appropriate at high temperatures above the 
Debye temperature (1076 K [Anderson, 1998]) and 
are entirely suitable for the Earth's mantle. We then 
used our results, combining them with data for other 
components and phases where appropriate, to cal- 
culate density and seismic velocity profiles along a 
geothermal gradient for different compositional mod- 
els. These profiles were then compared with profiles 
from a reference Earth model in order to test the com- 
positional models' fitness for the lower mantle. 

2. Method 

MD simulations were performed on superccUs of 
540-2500 atoms of orthorhombic {Pbnm) Mg-pv using 
the nonempirical V IB potential. VIB is a Gordon- 
Kim type model [Gordon and Kim, 1972], where 



the potential is obtained by overlapping ionic charge 
densities which are computed using the local den- 
sity approximation (LDA) [Hedin and Lundqvist 



1971]. The total energy is a sum of (1) the self- 
energy of each atom, (2) the long-range electrostatic 
energy computed using the Ewald method, and (3) 
the short-range interaction energy, i.e., the sum of 
the kinetic, short-range electrostatic, and exchange- 
correlation overlap energies from the LDA. Free oxy- 
gen ions are not stable and are stabilized in VIB by 
surrounding them with Watson spheres with charges 
equal in magnitude but opposite in sign to the ions, 
e.g., 2+ spheres for O'^^. Interactions are obtained 



3 



for overlapping ion pairs at different distances with 
different Watson spliere radii for eacli oxygen. The 
interactions are fit with a 23-parameter analytical ex- 
pression as functions of the interatomic distance r 
and Ui = Zi/Ri, where [7^, Zi, and Ri are the Wat- 
son sphere potential, charge, and radius for atom i. 
For each oxygen i, Ri is adjusted to minimize the to- 
tal energy at each time step. This gives an effective 
many-body potential. The oxygen ions respond to 
changes in their environment. For example, they are 
compressed at high pressures relative to low pressures. 
Previous work on Mg-pv using the related PIB model 
gave an equa tion of state i n excellent agreement with 
experiment \ Cohen, 1987 1. Also, work using VIB on 
MgO has shown that this m odel accurately predict s 
EOS and thermal properties [Inbar and Cohen, 1995|. 

Nominally charged ions give semiquantitative, but 
less accurate, results than desired for Mg-pv. This 
is due to a small degree of covalent bonding relative 
to ionic bonding as revealed by electronic structure 
calculations of cubic {Pmim) Mg-pv performed us- 
ing the first principles LAPW method \Cohen ei 
ai, 19891. These calculations show that while the 



Mg is nearly a perfectly spherical Mg^+ ion, there is 
some charge transfer from O to Si and a small co- 
valent bond charge (involving ^ 0.1 electrons) be- 
tween Si and O. By varying the ionic charges on Si 
{3.1+ to 3.55-h) and O (1.70- to 1.85-) to account for 
the covalency, but otherwise using the same methods 
as before, we compared the resulting zero tempera- 
ture pressure-volume data to the LAPW results of 
Stixrude and Cohen [1993]. The best agreement was 
found with Si'^ ''+ and 0^-*~. Using these charges, 
the resulting potential gives excellent agreement with 
the octahedral ro tation energetics obtained using the 
LA PW method \Stixrude and Cohen, 1993 ; Hemlei 



an 4 Cohen, 1996| (Figure g). 

MD runs were performed for 20 ps with a 1 fs 
time ste p using a six th-order Gear predictor-corrector 
scheme \ Cear, 1971 ] in the constant pressure-tempera- 
ture ensemble using the thermostat and barostat of 
Martyna et al. [1994]. Initial atomic positions for 
the genesis run were the same as those for the unro- 
tated octahedra of Stixrude and Cohen [1993]. Subse- 
quent runs at higher pressures or temperatures used 
the positions generated by a previous run. Statistical 
ensembles were obtained in ^2000 iterations for the 
genesis run and in < 1000 iterations for subsequent 
runs. 

We fitted pressure-volume-temperature (P-V-T) 
data obtained from the MD simulations to a thermal 



EOS of the form 

P{V, T) = Pa{V, To) + APth(r) 



(1) 



with a reference temperature Tq = 300 K and the 
thermal pressure APth relative to the 300 K isotherm. 
We analyzed our results using the third-order Birch- 



Murnaghan isothermal EOS [Birch, 1952] 



(2) 



with the Eulerian strain variable / and the coefficient 
f given by 



2/3 
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and the Universal EOS IVinet et al, 1987] 
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Pq = 3Kto {1 — x)x exp 



(X^„-l)(l-x) 



(3) 



, (4) 



where x — {V/VoY^^, Kt is the isothermal bulk mod- 
ulus, and K!p is its pressure derivative. 
APth is given by 



APth = 



dP 

To \dTyy 



dT= I {aKT)dT, (5) 

To 



where a is the volume coefficient of thermal expan- 
sion. If aKx is independent of T, then the thermal 
pressure can be represented by a linear relationsh ip 
\ Anderson, 1980, 1984 ; Anderson and Suzuki, 1983|] 



AP 



th 



a + bT, 



(6) 



which applies to a wide range of solids at high T, in- 
cluding minerals, alkali metals, and noble gas solids 



{Anderson, 1984 


]. The hnearity of APth 


at 300 K ] 


Masuda and Anderson, 1995 



als. There should be a small anharmonic correction 
at high T, which results in an additional cT^ term 
in (6), which is often sufficiently small so that it can 
be neglected [Anderson, 198^^]. Indeed, including the 
cT^ term did not statistically improve our fits. Like- 
wise, our fits were not improved by including volu- 
metric compression terms [e.g., Jackson and Rigden, 
1996]. 

Once the P-V-T data were fit, other parameters, 
such as a, 

p 



dT 
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its volume dependence, the Anderson- Griineisen pa- 
rameter, 

I din a 

Ot 



d\nV 



(8) 



the Griineisen parameter. 
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aKrV 



and its volume dependence, 

9 In 7 
dXnV 



(9) 



(10) 



were obtained by numerical differentiation. Cp was 
found by calculating enthalpy ( C/ + PV ) at each MD 
run point and differentiating with respect to T. Using 
the relation Cp = Cv(l + 7Q;T), equation (9) can be 
rewritten as 



( 



yaKrV 



Ta 



(11) 



3. Results 



-VlD calculaLiuus weie peifuiuied aL 40 piessuie- 
temperature (P-T) points, ranging from 0-140 GPa 



an d 300-3000 K, Lh 



of l!he luv 



le luwei iiiaiiLle ^^lauie >i.xiai iaiius uju aiii 
c/a follow smooth trends at lower mantle pressures 
(Figures |^a-3b), with deviations away from the Pbnm 
structure seen at P < 12 GPa and T > 2400 K. At 
3000 K and 12 GPa the structure can be observed 
to move towards cubic lattice parameters {b/a — > 
1, c/a — > \/2)- This point, in particular, is close 
to the melting curve of Mg-pv 
loz, 1987] ; ^nittle 



LIS cuveiiug 
(Table | l]). Axial 



the P-T cuudiL 



lUllS 

b/ti and 



I Heinz and Jean- 



d Jeanloz 198E; Poirier, 1989 1 



and such temperature-induced phase transitions have 



been predicted for Mg-pv by MD | Wolf and Bukowin- 
ski, 1987| and lattice dynamics [Matsui and Price 
1991[ . In addition, it has been found for RbCaPa and 



KMnFa perovskites by MD simulatio ns [Lewis ana 



Lepine, 198£ ; Nose and Klein, 1989 1 and observed 



experimenta lly in the fluoride perovskite neighborite 
(NaMgFa) ^Chao et al, 196l| . Regardless of these 



chai ngoD, the volumoa arc well behaved ao functionc of 



pre bouro and temperature (Figure 3c), and no other 
structural anomalies or melting was encountered dur- 
ing runs at other P-T points. 

The resulting P-V-T EOS fits and their reduced 
values are given in Table ||. Both the Birch- 
Murnaghan and Universal EOSs are in excellent agree- 
ment with experimental results. As for the accuracy 



of our thermal pressure expression, RMS differences 
between volumes found using equation (1) for either 
EOS and those found using isothermal EOSs are of 
the order of 10^'^ A'^. As for the differences between 
the two equation of state forms, Jeanloz [1988] com- 
pared them over moderate compressions and found 
that they are quite similar. More recent work by Co- 
hen et al. [2000] supports this, but they found that for 
large (>30%) strains and for determining parameters 
such as Vo, I^Ta^ and K'tq: the Universal EOS works 
best. So, for the rest of our analyses we use that here. 
However, given that under the pressure range studied, 
compressions will be no greater than 25%, so we can 
confidently compare our results with those of others 
that were found with the Birch-Murnaghan EOS. 

Once we determined the EOS parameters, we de- 
termined volumes for Mg-pv over a wide P-T range 
(-10-150 GPa, 0-3300 K). Using this extended data 
set, we determined the isothermal bulk modulus [Kt = 
— V{dP /dV)T] at those points and a, 5t, 7, and q 
using equations (7), (8), (11), and (10). We found 
{dK/dT)p=o = -0.0251 GPa K"!, cl ose to thos e 
fomid experimentally, -0.023 GPa [[VFawff et al 

1994|] and -0.027 GPa K^i 



[Fiquet et al, 1998 1, as 



well as one found by the inversion of multiple experi- 
mental data sets, -0.021 GPa K^^ [Jackson and Rig 



den, 1996 1 



Figure |j shows a at pressures from to 140 GPa. 
Experimental results shown have higher T slopes, as 
well as higher extrapolated values for most tempera- 
tures, though an average value found by Jackson and 
Rigden [1996] at zero GPa of 2.6 x 10"^ K'^ over 
300-1660 K matches our value over the same T range 
exactly. Kato et al. [1995] found an average value of 
2.0 ± 0.4 X 10^5 K-i over 298-1473 K and 25 GPa, 
close to our average value of 1.89 x 10~^ at the 
same conditions. The lack of quantum effects in our 
MD results can be seen at low temperatures; instead 
offending toward zero at zero temperature, a remains 
large. The 5t increases with T, but decreases as P 
increases (Figure ||) . The dependency of St on T de- 
creases as a function of P, with convergence to a value 
of 2.87 at ^^130 GPa, similar to the high-pressu re be- 
havior in MgO fo und by the same MD method [ Inbar 
and Cohen, 1995 1. Our value of St = 3.79 at ambient 
conditions is in excellent agreement with the experi- 
mentally derived value of Wang et al. [1994] of 3.73, 
though it is lower than the value of 4.5, found by 
Masuda and Anderson [1995] from the experimental 
data of Utsumi et al. [1995]. Other theoretical de- 
terminations of St are much higher, however. Hama 



and Suito [1998a] found a value of 5.21, on the basis 
of calculations using the LAPW results of Stixrude 
and Cohen [1993] and the Debye approximation for 
lattice vibration. Using MD and semiempirical po- 
tentials ^atsui, 1988 1, Belonoshko and Dubrovinsky 
[1996] found 5t = 5.80, and Patel et al. [1996] found 
5t — 7.0 using the same potentials, combining MD 
with lattice dynamics. Anderson [1998], using Debye 
theory, calculated zero pressure values ranging from 
4.98 at 400 K to ~4 at 1000 K < T < 1800 K, coming 
close to our values at those temperatures. 

The Griineisen parameter also increases as a func- 
tion of T and decreases as a function of P and ap- 
proaches a value of 1 at pressures of 130-140 GPa 
(Figure Our value of 1.33 at zero P and 300 K 
matches well with Wang et o/.'s [1994] value of 1.3 
and Utsumi et aZ.'s [1995] and Masuda and Ander- 
son's [1995] value of 1.45. Stixrude et al. [1992], us- 
ing the experimental data of Mao et al. [1991], found 
a higher value of 7 = 1.96. Values determined by 
the inversion of multiple experimentally determined 
P-V-T data sets match well also: 1.5 \Bina, 1995], 
1.33 \ Jackson and Rigden, 1996t , and 1.42 \ Shim ana 
Du jy, 2000]. Two shock wave studies find a value of 



1.5 [Duffy and Ahrens, 1993; Akins and Ahrens, 1999| ] 



on the basis of limited data sets: four for the former 
(with q fixed equal to 1) and two for the latter. Values 
of 7 from theoretical studies range from very close to 



ours, 1.279 ]ffama and Smto, 1998a and 1.44 {Hem- 


ley 
sky 


et al, 1987|, to 1.97 J 


Belonoshko and Dubrovin- 


1996; 


Patel et al., 1996]. Anderson's [1998] zero 



pressure, 300 K value of 1.52 is somewhat higher than 
ours, but at 400-1800 K his values of ^1.4 are very 
close to ours. 

The volume dependence of 7, q, is 1.03 at 300 K 
(Figure Our equation of state form, with APth 
linear in T and independent of volume, constrains q 
= 1 if the isochoric heat capacity Cv — SnR, the clas- 
sical harmonic value. As T increases, we find that q 
increases to a value of 1.07 at 3000 K owing to anhar- 
monicity. Values from inverted experimental data sets 
range from 1.0 \ Bina, 1995| to 2.0 \ Shim and Duffy, 
20(Q. Stixrude et al. [1992] found a high q value 
of 2.5 to go along with their value for 7. Akins and 
Ahrens'' s [1999] value was even higher, q = 4.0 ±1.0 
with Cv = 5ni?, but these are preliminary conclu- 
sions based on limited data. Patel et al. [1996] found 
a value of 3.0 at GPa, decreasing to 1.7 at 100 GPa 
using a combination of molecular and lattice dynam- 
ics. Anderson's [1998] Debye calculations gave val- 
ues close to unity at T > 1000 K, with q decreasing 



slightly with increasing T, to 0.82 at 1800 K. 

4. Discussion and Conclusions 

Taking sets of experi mental P-V-T data, we fit the 



high temperatur e data \ Fiquet et al., 1998, 2000 ; Sax- 



ena et al., 1999 and available experimental MgSiOa 
data to our thermal Universal EOS form (Table 2) in 
order to have consistent bases for comparison. The 
resulting fits have higher statistical uncertainties but 
compare well with our results. Volumes calculated 
along the lower mantle geothermal gradient of Brown 
and Shankland [1981], which has a starting T of 1873 
K at 670 km (Figure 1), show that those calculated 
from our EOS are 3 to 4 per unit cell larger, 
^--^2.5%, than those calculated from the inverted ex- 
perimental data sets (Figure ||). This corresponds 
to density differences of ~0.1 g cm~^. This is the 
opposite of the typical error of the local density ap- 
proximation to density functional theory, on which 
our method is based. The larger volume is due to 
the choice of ionic charges and other model assump- 
tions. If we fix Vo to the experimentally determined 
value of 162.49 [ |Mao et al., 199l| ] but otherwise 
use our EOS parameters, our calculated volumes are 
very close, within 1%, to the volumes derived from 
the inverted experimental data sets. The compar- 
isons also improve as pressure increases. To compare 
our model with other zero temperature theoretically 
derived EOSs, we calculated volumes at K from 0- 
140 GPa and fit an isothermal Universal EOS to them 
(Table |). 

Compositional models of the lower mantle can 
be tested against seismological models by examin- 
ing densities and seismic wave velocities as func- 



tions of depth. Candidates include pyrolite [Ring 



wood, 1975] and chondritic or pure perovskitc models. 



Studies support both the former [ Bina and Silver. 



199C; Stacey, 1996] and latter [Butler and Ander- 



son, 1978; Stixrude et al., 1992], though uncertainties 
in the thermodynamic parameters of the constituent 
minerals make it difficult to resolve this question with 
certainty. Indeed, several studies have been able 
to support both models depending on whether high 
(pure perovskitc) or low (pyrolite) values of a and 5t 
or 7 and q are used \Zhao and Anderson, 1994 ; An 
derson et al., 1995; Karki and Stixrude, 1999| [! 



Calculating the densities of Mg-pv along a geother- 
mal gradient [Brown and Shankland, 198l[ (Fig- 
ure ^), we see that they are ~0.2-0.3 g cm~'^ (^--^4- 
5%) too low compared with those from the Prelim- 
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inary Reference Earth Model (PREM) [ Dziewonsk-\ 
am', Anderson, 1981]. The calculated bulk sound ve- 
locities, — y^Ks/p, where Ks = Kt{1 + ^aT), 
are, in turn, ~0.4-0.6 km (^4-7%) too high corn- 



par ed~to~thF'akT35TsmsnwIogi^^ 



at, 1995; Montagner and Kennett, 1996| | (Figure 9b) 



Using the experimentally determined value of Vb = 
162.49 \Mao et al, 1991] in place of our calcu- 



lated value does increase the densities somewhat, but 
not enough to match the PREM values. 

As the perovskite found in the lower mantle should 
be a solid solution of the Mg and Fe end-members, 



we added iron by adjusting density by \Jeanloz ana 
Th impson, 1985(] , 



p(XFe) =p(0)(l+ 0.269 Xfc) 



(12) 



but we did not change the bulk moduli [ Mao e\ 
al. \ 1991] . Including 10 mol % Fe does cause den- 
sities to agree somewhat with the PREM values (Fig- 
ure p^), but the bulk sound velocities are still much 
too high (Figure 10b). Consequently, we find that 
the lower mantle cannot consist solely of (Mg,Fe)Si03 
o^okitc. 



peri8¥ 



'/Ve also Lesled pyroliLic composlLious cousisUug 
of mixtures of perovskite and magnesiowiistite (mw). 
Densities of mw were calculated using the thermody- 
namic data set of Fei et al. [1991], as well as Kts for 
the Mg end-member. Kss were then calculated using 
the relation 



a(P,T) =«(Po,T) 



V{Po,T) 



(13) 



with St and 7 from Inbar and Cohen [1995]. Iron was 
accounted for in Ks via the relationship \Jeanloz ana 
Th impson, 1983| ] 



Ks {Xf,) = Ks{0) + 17 Xf,. 



Xug = 0.93-0.96 for perovskite (pv) and 0.82-0.86 for 
mw \Katsura and Ito, 1996; Martinez et al., 1997|] 
were used. Densities and bulk sound velocities were 
calculated using Voigt-Reuss-Hill averaging for two 



siniphfied pyrolitc models (high [i] and low [i] Fe 
content) consisting of i vol % pv and 66 vol % mw 
and were found to be in excellent agreement with the 
seismological models (Figure |ll|). The two pyrolite 
models have partitioning coefficients -S^pc-Mg ^-^^ 
and 0.26, where 



T^pv— mw ^Fc/^Mg 



Fc— Mg vmw / vmw ' 
^Fo /^Mg 



(15) 



Experimental evidence suggests that for such a bulk 
composition, ^^pe-Mg' should increase in the mantle 
from ~0.20 (660 km) to -0.35 (1500 km), wi th Xp;^ 
decreasing and XJJJg increasing with depth \Mao e 



al., 1997]. Compositional models with high Fe con- 



tent of pv and low Fe content of mw (and vice versa) 
fall inbetween the results of those shown in Figure 11. 
Given the range of Fe-Mg partitioning between per- 
ovskite and magnesiowiistite under the appropriate 
pressure and temperature conditions, we find pyro- 
lite to be the most likely compositional model for the 
lower mantle. 

Given that other components (Al, Fe'^^) and phases 
(CaSiOs perovskite) should be present in the lower 
mantle, we do not expect exact agreement of a sim- 
plified pyrolite model with any seismological model. 
It is believed that, under lower mantle conditions, 
AI2O 3 is mainly incorporated i nto the Mg-pv struc- 
ture llrifune, l994l\Wood, 20"0^ . Generally, the effect 



of the incorporation of Al into Mg-pv on its physical 
properties has been considered small [e.g., Weidner 
and Wang, 1998]. However, Al, unlike Fe, causes 
significant distortion in the Mg-pv lattice ] 'Neil 



and Jeanloz, 1994], which may affect the compress- 



ibility. Recent experimental work by Zhang and Wei- 
dner [1999], at pressure of up to 10 GPa and temper- 
atures of up to 1073 K, indicates that compared with 
MgSiOs, Mg-pv with 5 mol % Al has a smaller Kt 
value (232-236 GPa) and a {dKT/dT)p value more 
than double that in magnitude. The values for ao; 
{da/dT)p, and 6t are also larger. Smaller bulk mod- 
uh and the higher density (0.2% at 300 K and GPa) 
would cause seismic velocities to decrease. In addi- 
tion, the presence of Al tends to equalize the par- 
titioning of Fe into perovskite and magnesiowiistite 
and may allow garnet to coexist with perovskite in 
the uppermost —100 km of the lower mantle ] Wooc 



and Rubie, 1996; Wood, 200C]. These partitioning ex- 
periments, however, were done at high /02, so it is 
uncertain how applicable they are to the mantle. The 
effect of the presence of Fe'^"'' in Mg-pv on the EOS 
and elasticity, while known for defect concentrations 
and electrostatic charge balance \McCammon, 1997] 



199S], is unknown. As for CaSiOs perovskite (Ca-pv), 



experimental data suggest that its elastic properties 
are in excellent agreement with seismological mod- 
els and thus would be invisible in the lower mantle 
[Wang et al, 199(;; Hama and Suito, 1998b] . 



Performing MD calculations over the range pres- 
sures and temperatures found in the Earth's mantle, 
we find a thermal EOS that is in excellent agreement 



7 



with experimental results. Thermodynamic parame- 
ters can be derived that agree well with experimen- 
tally determined values and that can be confidently 
interpolated to conditions found in the lower mantle. 
Moreover, no solid-solid phase transitions or melting 
was found during the runs under lower mantle condi- 
tions, so orthorhombic MgSiOa perovskite should be 
stable to the core-mantle boundary. Using these re- 
sults, we find that pyrolite with K^^Zug — 0.26-0.34 
is the most likely compositional model for the lower 
mantle. 
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Figure 1. Experimental pressure-temperature coverage of (Mg,Fe)Si03 perovskite. Open symbols are for static 
experiments on XMg = 1-0, solid symbols are for XMg = 0.9, and open symbols with crosses are shock wave data. 
Hexagons, Kudoh et al, [1987]; right-facing triangles, Ross and Hazen, [1990]; left-facing triangles, Wang et al, 
[1994]; upward pointing triangles, Fiquet et al., [1998]; downward pointing triangles, Fiquet et al, [2000]; diamonds, 
Saxena et al., [1999]; open circles, Utsumi et al, [1995]; solid circles, Knittle and Jeanloz, [1987]; squares, Funamori 
et al., [1996]; stars, Mao et al., [1991]; pentagons, Kato et al., [1995]; circle with cross, Duffy and Ahrens, [1993]. 
Solid curve is the lower mantle geothermal gradient of Brown and Shankland [1981]. Static experimental work has 
direct coverage of approximately the upper one third. The shock wave data point that falls near the geotherm is 
perovskite plus magnesiowiistite. 
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Fig. 2 
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Figure 2. Energetics of octahedral rotations: calculated total energies relative to the cubic structure {Pm2>m) 
using VIB (curves) and LAPW [Stixrude and Cohen, 1993; Hemley and Cohen, 1996] (symbols) as a function of 
(left) R and (right) M point rotation angles as represented by the fractional change in the oxygen coordinate {Sr 
and Sm)- The orthorhombic structure {Pbnm) occurs at 5r = 0.0912 and 5m = 0.0766 at P = GPa. 
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Fig. 3a 
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Fig. 3c 
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Figure 3. Structural parameters and volumes of MD results: (a) b/a axis ratio and (b) c/a axis ratio. At low 
pressures and high temperatures the structure begins to deviate from orthorhombic Pbnm, but no phase transitions, 
with the possibility of the structure heading toward cubic at 12 GPa and 3000 K. (c) Volume. Symbols are MD 
results and curves are Universal EOS fits every 300 K. 




Figure 4. Volume coefficient of thermal expansion. Solid curves are from this work. Other symbols, from 
experimental work as follows: dashed curve (0 GPa), Wang et al. [1994]; circle (0 GPa), Fiquet et al. [1998]; square 
(0 GPa) and dashed-dotted line (25 GFa),Funamori et al. [1996] 



16 




T(K) 

Figure 5. Anderson- Griineisen parameter as a function of temperature. The circle is the experimentally determined 
value from Wang et al. [1994]. The square is the value determined by Masuda and Anderson [1995] from the 
experimental data of Utsumi et al. [1995] . Open circles are 6ts determined from Debye theory by Anderson [1998] . 
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Figure 6. Griineisen parameter as a function of temperature. The solid circle is the experimentally determined 
value from Wang et al. [1994], and the square is the value determined by Masuda and Anderson [1995] from the 
experimental data of Utsumi et al. [1995]. Other symbols arc from inversion of multiple experimentally determined 
data sets: triangle, Bina [1995]; inverted triangle, Jackson and Rigden [1996]); and diamond. Shim and Duffy 2000]. 
Open circles are 7s determined from Debye theory by Anderson [1998] . 




Figure 7. The In 7 versus In V. The slope of each curve gives q: soUd curve and squares, 300 K; dashed curve and 
circles, 1200 K; dotted curve and triangles, 2100 K; dash-dotted curve and inverted triangles, 3000 K. The lower 
temperatures are all close to 1, with a slight rise to 1.07 at 3000 K. 
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Figure 8. Volumes of MgSiOa perovskite along the geothermal gradient of Brown and Shankland [1981]. Dash- 
dotted curve is volume from our Universal EOS fit. Dashed curve is from our Universal fit with Vq fixed to 162.49 

[Mao et aL, 1991]. Dotted curve is for experimental data sets with Vq fixed to 162.49 A"^ (sec Tabic 2). Both 
high temperature (F98-|-F00-|-S99) and all MgSiOa data sets fall on the same line. Shaded areas indicate the 
uncertainties in the fits of the experimental data sets. 



20 



Fig. 9a 
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Figure 9. (a) Densities and (b) bulk sound velocities of MgSiOa perovskite calculated along the lower mantle 
geotherm of Brown and Shankland [1981] using the thermal Universal EOS. Solid curve is for uncorrected volumes 
and dashed curve is for Vq set to the experimental value of Mao et al., [1991]. Squares are densities from PREM 
[Dziewonski and Anderson, 1981]. Circles are velocities from the akl35-f seismological model [Kennett et al, 1995; 
Montagner and Kennett, 1996]. 
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Fig. 10a 
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Figure 10. (a) Densities and (b) bulk sound velocities of (Mg,Fe)Si03 perovskite calculated along the geotherm 
of Brown and Shankland [1981] using the experimental value of Vq Mao et al, [1991] for = 0.00 (solid curve), 
0.05 (dashed curve), and 0.10 (dotted curve). Squares arc densities from PREM [Dziewonski and Anderson, 1981]. 
Circles are velocities from the akl35-f seismological model [Kennett et al, 1995; Montagner and Kennett, 1996]. 
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Fig. 11a 
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Figure 11. (a) Densities and (b) bulk sound velocities of two simplified pyrolite models (67 vol % pv, 33 vol % 
mw) along the geotherm of Brown and Shankland [1981]. Solid curve is = 0.93, XjJJg = 0.82. Dashed curve 
is = 0.96, XJJJ"^ = 0.86. Squares arc densities from PREM [Dziewonski and Anderson, 1981]. Circles are 
velocities from the akl35-f seismological model [Kennett et ai, 1995; Montagner and Kennett, 1996]. 
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Table 1. MD P-V-T Data Used in EOS Fits 



GlPa 


T K 


a, A 


b, A 


c, A 


V A3 





300 


4.8166 


A AO O O 

4.9283 


6.9649 


165.33 





600 


4.8299 


4.9397 


6.9787 


166.50 





900 


4.8439 


4.9519 


6.9944 


167.77 





1200 


4.8586 


4.9648 


7.0107 


169.11 





1500 


4.8732 


4.9769 


7.0289 


170.48 





1800 


4.8886 


4.9901 


7.0455 


171.87 





2100 


4.9109 


5.0145 


7.0578 


173.80 





2400 


4.9311 


5.0354 


7.0688 


175.52 





2700 


4.9525 


5.0538 


7.0831 


177.28 





3000 


4.9672 


5.0633 


7.1076 


178.76 


12 


900 


4.7790 


4.8856 


6.8895 


160.86 


12 


1200 


4.7917 


4.8951 


6.9024 


161.90 


12 


1500 


4.8036 


4.9051 


6.9165 


162.97 


12 


1800 


4.8169 


4.9157 


6.9288 


164.06 


12 


2100 


4.8295 


4.9256 


6.9441 


165.18 


12 


2400 


4.8429 


4.9367 


6.9586 


166.36 


12 


2700 


4.8628 


4.9539 


6.9690 


167.88 


12 


3000 


4.9035 


4.9531 


6.9695 


169.27 


25 


300 


4.6971 


4.8093 


6.7728 


152.99 


25 


600 


4.7073 


4.8169 


6.7830 


153.80 


25 


900 


4.7181 


4.8250 


6.7948 


154.68 


25 


1200 


4.7285 


A o o o r\ 

4.8330 


6.8055 


155.53 


25 


1500 


4.7390 


A A f\ r\ 

4.8409 


6.8166 


156.38 


25 


1800 


4.7497 


4.8487 


6.8291 


157.27 


25 


2100 


4.7609 


4.8572 


6.8408 


158.19 


25 


2400 


4.7728 


4.8658 


6.8514 


159.11 


25 


2700 


4.7850 


4.8736 


6.8645 


160.08 


25 


3000 


4.8015 


4.8859 


6.8773 


161.34 


50 


900 


4.6194 


4.7307 


6.6445 


145.20 


50 


1200 


4.6280 


4.7366 


6.6531 


145.84 


50 


1500 


4.6366 


4.7414 


6.6632 


146.48 


50 


1800 


4.6450 


A n A n A 

4.7474 


6.6726 


147.14 


50 


2100 


4.6545 


4.7524 


6.6821 


147.81 


50 


2400 


4.6627 


4.7589 


6.6920 


148.49 


50 


2700 


4.6728 


4.7644 


6.7014 


149.20 


50 


3000 


4.6814 


4.7700 


6.7129 


149.90 


140 


300 


4.3572 


4.5023 


6.2698 


123.00 


140 


600 


4.3627 


4.5046 


6.2768 


123.35 


140 


900 


4.3687 


4.5076 


6.2839 


123.74 


140 


1200 


4.3747 


4.5097 


6.2903 


124.10 


140 


1500 


4.3806 


4.5121 


6.2967 


124.46 


140 


1800 


4.3867 


4.5146 


6.3027 


124.82 


140 


2100 


4.3923 


4.5164 


6.3104 


125.18 


140 


2400 


4.3986 


4.5186 


6.3170 


125.55 
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Table 1. (continued) 



P, GPa 


T, K 


a, A 


b, A 


c, A 


V, A3 


140 


2700 


4.4048 


4.5209 


6.3235 


125.92 


140 


3000 


4.4109 


4.5230 


6.3309 


126.31 
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Table 2. Equation of State Parameters'* 



Vo, A3 KTo,GPa K^^ a, GPa b, GPa K"! x^ GPa 



This work 



Birch-Murnaghan EOS 


165.40(6) 


274(1) 


3.73(2) 


-2.00(2) 


0.00667(5) 


0.08561 


Universal EOS 


165.40(5) 


273(1) 


3.86(2) 


-1.99(1) 


0.00664(5) 


0.06712 


Experimental results 














M91 


162.49(7) 


261(4) 


4b 








F98 


162.65(15) 


262(6) 


3.41(22) 


-1.79(21) 


0.00597(70) 


1.47438 




162.49^= 


267(3) 


3.25(10) 


-1.79(21) 


0.00597(70) 


1.41566 




162.60(9) 


264<l 


3.33(2) 


-1.79(21) 


0.00596(69) 


1.41333 




162.49^= 


264^^ 


3.41(3) 


-1.80 35 


0.00601(73) 


1.36143 


FOO 


162.74(56) 


239(8) 


4.41(18) 


-1.85(6) 


0.00618(20) 


0.61052 




162.49^^ 


245(4) 


4.29(11) 


-1.85(6) 


0.00618(20) 


0.59369 




161.68(18) 


264^ 


3.88(2) 


-1.85(6) 


0.00617(20) 


0.60570 




162.49^= 


264^ 


3.52(4) 


-1.73(8) 


0.00577(28) 


1.23307 


S99 


163.25(63) 


242(6) 


4.45(10) 


-1.72(1) 


0.00574(2) 


0.26719 




162.49^= 


256(5) 


4.20(10) 


-1.72(1) 


0.00572(2) 


0.26227 




161.10(28) 


264^ 


4.06(4) 


-1.71(1) 


0.00572(2) 


0.26645 




162.49^^ 


264<l 


3.89(2) 


-1.75(2) 


0.00584(7) 


0.31819 


F98 + FOO + S99 


162.93(33) 


235(6) 


4.72(13) 


-1.85(12) 


0.00618(38) 


1.56307 




162.49^= 


244(4) 


4.52(13) 


-1.86(12) 


0.00621(40) 


1.55478 




161.69(22) 


264'! 


4.08(1) 


-1.87(12) 


0.00622(38) 


1.61901 




162.49^= 


264^^ 


3.75(2) 


-1.73(15) 


0.00577(7) 


2.22349 


All 


162.31(10) 


249(6) 


4.40(20) 


-1.92(10) 


0.00639(32) 


0.60337 




162.49^^ 


243(8) 


4.55(28) 


-1.90(10) 


0.00635(34) 


0.61628 




162.01(22) 


264^ 


3.92(8) 


-1.88(10) 


0.00628(38) 


0.67860 




162.49^= 


264^ 


3.75(2) 


-1.67(20) 


0.00557(65) 


1.04352 



''MD are best fit EOS parameters to molecular dynamics P-V-T data. Experimental results are the best fit Universal 
EOS parameters to experimental data sets, except data taken from Mao et al. [1991] (M91). Remaining experimental 
results are F98, Piquet et al. [1998] (27 data points); FOO, Piquet et al. [2000] (38 data points); S99, Saxena et al. [1999], 
(37 data points); All, all XMg = 1.0 data (363 data points; see Figure 1). Values in parentheses are standard deviations. 

''Fixed equal to 4. 

'^Fixed to Uo from Mao et al. [1991]. 
''Fixed to Ktq from Yeganeh-Haen [1994]. 
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Table 3. Zero Temperature Equation of State Parameters 



V, A3 K, GPa K' 



This work* 164.22 280 3.84 

PIB'^ 164.78 252 4.05 

LAPW^ 160.74 266 4.2 

PWPPd 157.50 259 3.9 



^Universal fit. Other theoretical fits are third-order Birch- 

Murnaghan fits. 

Cohen [1987]. 
'^Stixrude and Cohen [1993]. 
"^Wentzcovitch et al. [1995]. 



